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SUMMARY 

Flux-vector and flux-difference splittings for the inviscid terms of the compressible flow 
equations are derived under the assumption of a general equation of state for a real gas 
in equilibrium. No unnecessary assumptions, approximations or auxiliary quantities are 
introduced. The formulas derived include several particular cases known for ideal gases and 
readily apply to curvilinear coordinates. Applications of the formulas in a TVD algorithm 
to one-dimensional shock-tube and nozzle problems show their quality and robustness. 
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1. INTRODUCTION 


Several split-flux formulas are known for upwind differencing of the inviscid terms in 
CFD codes for ideal gases. The generic formula is due to Gudonov [l] and is based on the 
exact solution of Riemann’s initial-value problem, representing the interaction of two fluid 
parcels by finite-amplitude waves. 

Numerical efficiency justifies the introduction of approximations to the Riemann solu- 
tion, which leads to various simplifications of the flux formula, accompanied by considerable 
savings in computational expenses. The most popular “approximate Riemann solvers” are 
the flux-vector splittings by Steger and Warming [2] and by Van Leer [3] and the flux- 
difference splittings by Roe [4] and by Osher [5]; see e.g., the review by Harten, Lax and 
Van Leer [6]. 

With the current interest in high-temperature flows, real-gas effects* must be included, 
requiring appropriate modifications of all of the above split-flux formulas. Colella and Glaz 
[7j extended the numerical procedure for obtaining the exact Riemann solution to the real- 
gas case. Grossman and Walters [8], as well as Vinokur and Liu [9] extended the formulas 
of Steger and Warming, Van Leer and Roe, while Glaister [10] presented an extension of 
the Roe splitting. 

As regards the approximate Riemann solvers for real gases, the derivations in [8] 
are the least general, introducing unnecessary assumptions, approximations and auxiliary 
quantities. A more careful and comprehensive analysis is presented in [9]. Yet the combined 
formulas in [8-10] far from exhaust the possibilities, this forms the chief motivation for the 
present paper. 

The present derivation of split-flux formulas has the following features: (1) it includes 
several particular formulas derived elsewhere for ideal gases, (2) it avoids unnecessary 
assumptions or approximations, (3) it avoids unnecessary auxiliary quantities, and (4) it 
readily extends to curvilinear coordinates. 

In the next section we briefly discuss the equation of state (EOS) for the real gas 
and some related thermodynamics quantities. The detailed derivation of split-flux for- 
mulas is given in Section 3. Application of this construction in a TVD algorithm [ll] is 


* In this paper the definition of “real gas” is broader than that conventionally used in 
thermodynamics, by referring to any gas that is not both thermally and calorically perfect. 
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demonstrated for one-dimensional shock tube and nozzle problems. 


2. EQUATION OF STATE 

We begin by assuming that the real gas is described by the general equation of state, 

P = t = 1,2,3... (2.1) 

where p, p, e, Y, are, respectively, the pressure, density, and specific internal energy of 
the gas and the mass fraction of the i-th species in the gas. In this paper, we will restrict 
ourselves to the gas in chemical equilibrium, i.e., 

P = p(p,e), (2.2) 

so that derivatives with respect to Y, disappear, simplifying the algebra. Yet, the con- 
cepts and algebraic steps needed to describe a non-equilibrium gas are very nearly the 
same as for the equilibrium gas. In fact, some of the special care in formulating the EOS 
of an equilibrium gas can be relaxed when including nonequilibrium (finite-rate reaction) 
effects, because the overall numerical procedure becomes more straightforward. The com- 
putational effort, of course, is increased enormously. 

For an ideal gas, ( 2 . 2 ) reduces to p = (7 — l)pe, where 7 is the ratio of enthalpy h to 
internal energy e, and is a constant. 

The speed of sound is 

a 2 = pp + ppjp 2 . (2.3) 

Here and in this paper, p p and p e denote the partial derivatives of p with respect to p and 
e while holding other variables fixed. 

Several authors [7-9] have advocated the use of an equivalent 7 for real gases. As in 
the ideal gas, we define 

7 = h/e, h — specific enthalpy; (2.4) 

combining this with the definition, h = e + p/ p, gives 

P = (7 - 1 )pe. (2.5a) 

This is identical in form to the EOS of the ideal gas, but now 

7 = Tf(p»e)- (2-56) 
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We note that, although the range of variation of 7 may be limited, it may be of great 
importance. Since it is the factor (7 — 1) that appears in the equations, a small change in 
7 can result in an appreciable change in p, e.g., 7 = 1.2 doubles the pressure from 7 = 1.1 
even though A7 = 0.1. On the other hand, changes in flow variables of several orders of 
magnitude can occur with 7 fixed, as in the case of an ideal gas. Note that for a real gas 
7 in general will not be a monotone function of temperature. 

The speed of sound expressed in terms of 7 now has the form 

a 2 = IP /p + 1 [(7 - l)e7e + P7p] , (2.6a) 

where 7 e = d^/dt and 7 P = d^/dp. Since for an ideal gas (2.6a) reduces to a 2 = 7p/p, 
this motivates the introduction of an auxiliary variable T such that 

a 2 = Tp/p = T(7 - l)e, (2.66) 


where T = T(p, e). 

We caution that the appearance of these different “equivalent 7 V* may at times add 
more confusion than insight. 

From thermodynamic principles, it is possible to calculate any thermodynamic variable 
for each pair of state quantities (p, e). In practice, a table or a least-square fitted surface 
is generated a priori for reasons of computational efficiency. We adopt the latter approach 
because interpolation is avoided, resulting in greater efficiency. However, care must be 
taken to ensure, while fitting, not only that the error is kept whithin reasonable bounds 
but also that no numerical oscillations are introduced. 


3. CONSTRUCTION OF SPLIT FLUXES 


To illustrate how the three split-flux formulas are constructed, we first consider the 
ID Euler equations, 

where U* = [p,pu,pE], and F* = [pu,pu 2 + p, ( pE + p)u] . 

As indicated previously, the EOS now is expressed in terms of all variables U, viz. as 
p = p(p(U), e(U)). Thus the Jacobian matrix is readily derived, 


A = 


dF 

au 


— A e 4- A p, 


(3.2a) 
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where 


A e 


/ o i o A 

-u 2 {2- p e /p)/2 + {p-p e e)/p u(2 — p e /p) Pe/p 
\u\-H + U 2 p e /2p + {p-p e e)/p] H-u 2 p e /p u(l +p e /p)J 


contains the derivative p e , and the matrix containing p p is 

/ 0 0 (A 

A p = Pp-p/p 0 0 1, 

V u (Pp - p/p) 0 °J 

where H — E + pjp is the specific total enthalpy. 

The eigenvalues of these matrices are respectively: 


(3.26) 


(3.2c) 


A(A) = u — a, u, u + a. 

(3.3a) 

A(A e ) = u - a e , u, u + a e . 

(3.36) 

A(Ap) = 0, 0, 0, 

(3.3c) 

where a 2 = p(p e /p + l)/p = a 2 — (p p -p/p). Thus the matrices A and A e have a complete 
set of eigenvectors, but A p does not. We note that a e = a and A e = A for a gas in which 
p depends linearly on p, viz. a thermally perfect gas. Furthermore we find 

F = F* + F\ 

(3.4) 

where 

F h = AU, 


and 

F' - A P U. 


That is, the flux vector F no longer posseses the property of homogeneity, rather is a 
sum of homogeneous and inhomogeneous parts. Since the matrix A has a complete set of 


eigenvectors, it can be readily diagonalized by a similarity matrix S whose column vectors 
are the right eigenvectors of A: 


A = SAS-‘, diagA = A(A). 


(3.5) 


§3.1 Steger- Warming Splitting 
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By splitting the eigenvalues, 


A = A + + A - , (3.6) 

the real-gas version of Steger- Warming flux-vector splitting is obtained with F^ = F^ +Fj^. 
The split fluxes are suited for upwind differencing. Note that, unlike in [8], the true speed 
of sound a of the gas is used to determine the switching in (3.6). Since A p does not have a 
complete set of eigenvectors, an equation that includes F alone does not have a hyperbolic 
charater. In consequence, central differencing, corresponding to the pseudo splitting 

F ± =F ± fc +^F, (3.66) 

may be appropriate. This has been verified in the numerical experiments of Section 4. 

An entirely different approach, based on the “beam scheme” [12], is described in [9]. A 
discussion and numerical comparison of the various real-gas versions of the Steger- Warming 
splitting is of limited interest, since this formula, among all split-flux formulas, has been 
shown in various studies (see e.g. [13]) to yield the least accurate solutions in the ideal-gas 
case. The formula derived above is as efficient as any real-gas version can be and appears 
satisfactory in practice. 

One cause of inaccuracy in the Steger- Warming splitting is that the split-flux deriva- 
tives are discontinuous when any of the eigenvalues A (A) changes sign. This gives rise to 
“glitches” in numerical solution, unless some extra smoothing is built into the flux for- 
mula. To solve this problem radically, Van Leer [3] developed a continuously differentiable 
splitting. 


§3.2 Van Leer Splitting 


In this section we present a most general and compact derivation of the Van Leer 
splitting that includes a family of flux choices and is easily applicable to curvilinear co- 
ordinates. It further is independent of the EOS used, and does not require homogeneity 
of the flux vector. Splitting takes place only when there exist eigenvalues of mixed signs, 
i.e., as M 2 < 1 for the system (3.1) with eigenvalues (3.3a); M is the local Mach number. 
Since u = Ma, the mass flux F j — paM is split as: 

F l =F++Fr, 

F± = ±— pa[M ± l) 2 . 

4 

Note that the split fluxes have vanishing slope as M — * ^1, yielding smooth switching. 
This forms the basis for splitting of the remaining fluxes, namely, for expressing F 2 and 
F 3 in terms of F±. Let us write 


F 2 = pu 2 + p 
= ( pua)M + p , 


(3.8a) 
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where the first term on the RHS, representing the convection of momentum, already has 
the desired factor, 

M = [(M+ l) 2 - (M - l) 2 ]/4. (3.86) 

Now we must see if p can be written in terms of (M + l) 2 and (M — l) 2 . Let us try the 
combination, 


p = x{(M + l) 2 - (M - l) 2 ] + y[(M + l) 2 + (M - l) 2 ] 
= 4xM + 2 y(M 2 + 1). 

This leaves us one equation for two unknowns x and y. By choosing 

2y = p, 


which must be true anyway for M = 0, we find 

x = — -pM 

4 


for arbitrary M. Thus p is recast as 

p = -^[F+(-u + 2a) + Ff(-u - 2a)]. (3.9) 

pa* 

Substitution in (3.8a) yields the splitting of Fz, 

= f 2 + + Ff, 

Ft =F, ± (u- ~(tiT2a)] (3.10) 

r'* 

= F±[u - (u T 2a)/r]. 


We note that the last equality is obtained by the use of definition (2.66) and is identical 
to the formula for an ideal gas. 

We turn now to the splitting of the energy flux F 3 . Again the flux contains a convection 
term: 

Fz — puE + pu 
= ( paE)M + pu, 

where the factor M again can be represented as in (3.86). Since p has factors (u ± 2a), it 
is natural to assume a quadratic function in u for pu; for symmetry reason this must have 
the form, 

pu = Fj + (/u 2 + 2mua + na 2 ) + Ff ( lu 2 — 2mua + na 2 ) 

= (/ + m)pu 3 + (m + n)pua 2 . 
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This leaves us many possibilities. We choose to eliminate pu 3 by letting 


/ + m = 0 , 


hence 


m + n = 


P 



Consequently, a family of oo choices, with single parameter m, results for F3. 


Fz - + F 3 " , 

F 3 ± =F 1 ± [ff-m(uTfl) 2 l ) 


(3.12) 


The identical result is obtained when letting (m + n) vanish, instead of (/ + m), because 
of the symmetric occurrence of u 2 and a 2 . Van Leer’s splitting is a member of this family, 
found by requiring that the terms in the square bracket form a perfect square. This leads 

to 

h 


m = 


a 2 

1+2& 


(3.13a) 


and 


1 ± (ua ± 2h) 

V 1 (a 2 + 2 h) ' 


(3.136) 


For an ideal gas, this results in the nice property that one eigenvalue of each split-flux 
Jacobian vanishes, because F± depends solely on F± and F±: 


F* = const 



(3.13c) 


As a result, numerical diffusion is minimal [14] and sharp steady shocks can be obtained. 
This property, however, is lost for the real gas. The particular eigenvalues still are very 
close to zero and may be either positive or negative. In the former case the formula has 
slightly increased dissipation, in the latter cause it may actually make upwind differencing 
unstable. Neither effect is noticeable in practice; see numerical tests in Section 4. Another 
obvious choice is 


m = 0, 

= FfH. 


(3.14) 


This has the advantage of being the simplest and most efficient formula possible; the ideal- 
gas version has one eigenvalue close to zero and positive. It was derived independently by 
Hanel et al [15]. 
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In the derivation of Vinokur and Liu [9], the one-parameter family of formulas is 
governed by the parameter xp in the split energy flux: 


Ft = f ' ± { l(r 2 ( p , ‘- ± i) 2lll ~ + < A " r^T) + *(» t “) 2 }• ( 315 “) 


The first term corresponds to Van Leer’s flux for an ideal gas, the second term vanishes 
for an ideal gas. Vinokur and Liu [9] therefore choose r/j = 0 as the extension of Van Leer’s 
scheme to the real-gas case. Our extension (3.136) corresponds to 

(r + i)(i + 2£)' 

while (3.14) is generated by 

* = (3-15c) 

Clearly, the two formulas (3.12) and (3.15a) generate the same split energy fluxes. Our 
formula, however, does not favor the original Van Leer flux and therefore is more convenient 
in searching for other special cases. 

When extending the split fluxes (3.10) and (3.12) to curvilinear multidimensional 
moving coordinates, we benefit from the fact that p and pu have been expressed in terms 
of F^. The derivation is straightforward and yields the following formulas for the split 
fluxes in the direction of a coordinate (: 


(3.156) 


= ±--(u±a) 2 = ±\pa{M± l) 2 , 

4 a 4 

Ft = F?\u - (,(u T 2a)/r], 
f? = F 1 ± |«-f,(«=F2a)/r], 

Ft = Ftfw - (,(u T 2o)/r], 

Ft = Ft\H + (,(u T 2a)/r - m(i =F 2a) 2 ]. 


Here u, v and w are the Cartesian velocity components and u is the contravariant 
^-velocity, i.e. the velocity normal to a moving surface of constant £; M is the correspond- 
ing contravariant Mach number. The form is the same as for the ideal gas except the T, 
instead of 7 , is used (see Thomas, Van Leer and Walters [16]). 

Note that none of the differentiable splittings derived above involves special derivatives 
of the EOS. The only thermodynamic quantities needed are p and a, which are standard. 
Unlike in [ 8 ], no assumption of homogeneity of the flux F is required. 
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§3.3 Roe Splitting 

A 

To construct Roe’s flux-difference splitting, one usually defines an average state U 
such that 


AF = AAU, (3.17a) 

A = A(U), (3.176) 

U = U(Ui,,Ur), (3.17c) 


where 


and the matrix A, 


A(.) = (.)k-(*)l, 


A(U) = ^ 


= I — u 2 (2 - Pe/p)/ 2 -I- (Pp - p e e/p) 
u[-H + u 2 p e /2p + (pp - p e e/p)] 


u(2 -p e /p) 
if - u 2 p e /p 



(3.18) 


has precisely the same form as in (3.2a)-(3.2c). Note that the difference between states 
“L” and “R” is not necessarily small. 

The object now is to find an average state such that (3.17a) is satisfied exactly for all 
admissible pairs (U£,,Ufl). In the ideal-gas case this is easily accomplished since (3.17a) 
represents three relatively simple equations with three unknowns, the components of U. 
(Note, for instance, that the term (pp -p e e/p) appearing in the first column of A vanishes.) 
In fact, the matrix A(U) is completely determined by only two average state quantities as 
the density does not occur explicitly. The average Jacobian depends only on the quantities 
u and H\ an average density p may be chosen freely, although there is an obvious choice 
for it (see (3.20c)). 

For a real gas, the simplicity is lost. The nonlinear system (3.17a) must be solved 
numerically and in general may have more than one solution or no solution at all. In order 
to find a practical formula for A we must relax the constraint (3.176) and allow independent 
averages of more than three state quantities to enter the elements of A. Judging from the 
form of A in (3.18) we have six non-constant elements to play with, so there is a room 
for six independent averages. We shall choose the set ( u,p,e,H,p e ,p p ). We note that, 
the relation of the total enthalpy to the other flow quantities, i.e. H = e + p/p + u 2 / 2 
while holding pointwise, does not necessarily hold for the average state quantities which 
are functions of two states. Similarly, while p = p(p,e), p is not obtained by substituting 
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{p, e ) in the EOS, but rather must be defined separately; this is also the case for p p and 
p e . In the following derivation there is no need for an explicit definition of p. 

The mass-flux equation in (3.17a) is automatically satisfied for any average state. To 
satisfy the momentum-flux equation, we choose 

(i) A pu — pAu + uA p; (3.19a) 

(ii) A pu 2 = 2puAu + u 2 A p. (3.196) 

Let the Roe-average operator p be defined as: 

M)= TRl * + / Lh , r = pi. ' (3.20a) 

fR +r L 

The conditions (i) and (ii) are met by 



a = p(u), 

(3.206) 


P = r L r R . 

(3.20c) 

Assuming further that, 

(iii) Ape = pAe + eAp, 

(3.21a) 

which is met by 

e = p(e), 

(3.216) 

the second equation of (3.17a) becomes 



RHS 2 = A pu 2 + (p e Ae + p p Ap) 

LHS 2 = A pu 2 + A p. 

Clearly, the condition to be satisfied is 

(iv) A p = p e Ac + p p Ap. (3.22) 

Finally, we find that the energy-flux equation is readily satisfied by setting 

(v) A puH = pH Au + uA pH, (3.23a) 

hence 

H = p{H). (3.236) 
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Left to be completed are the definitions of p e and p p . Glaister [10] suggests 


Pe = 

Pp = 


Pe ~ Pw 
Ae 

Ps ~ Pi v 
Ap 


1 


(3.24a) 


where Pe,Pw ,Ps, and p ^ are approximated at points shown in the sketch by algebraic 
averaging 



i.e., 


Pe = ~(pr + Pb), Pb = p(PL,e R ), 
Pw = ~[PL+Pa), PA~P{PR, e L), 
Pn = 2 (pl + Pb), 

Ps = ^(Pfl + Pa)- 


(3.246) 


With this choice, (3.22) is met precisely. However, the artificially introduced states, A and 
B, cause problems at discontinuities, in particular for a contact discontinuity across which 
the density and internal energy jump appreciably. If the EOS is a non-convex function(as 
for an ideal gas) and L and R are two admissible discrete states, then A and B can lie 
outside the region of validity of the EOS. This will manifest itself in a calculation in the 
form of oscillations near the discontinuity, and/or an inaccurate solution. 


Here we propose another formula that uses derivative information at the average 
state (p, e), which is guaranteed to lie between states L and R by virtue of (3.20)-(3.2l). 
Introducing 

Pp — Pp{p, e), Pc=Pe{p,e), 
we divide the residual of (3.22), 


6p = Ap - PpAp - p e Ae, 
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equally over the terms, i.e., 


PpAp = PpAp + 0.56p, p e Ae = p e Ae + 0.5 6p. 


This leads to 

Pe = r[(Ap - PpAp)/ Ae + p e ], 

\ (3-25) 

Pp = 2 ((^P “ PeAe)/Ap + Pp). 

It is easily verified that, of all the pairs ( p p ,p e ) satisfying (3.22), the pair (3.25) lies closest 
to (p)>, pie) and must be considered the optimal choice regarding consistency. The formulas 
suggested by Vinokur and Liu [9] are based on a different choice of independently averaged 
quantities, with different constraints, and are more complicated. 


4. NUMERICAL TEST 

Extensive tests over a wide range of flow conditions have been conducted to validate 
the accuracy of the present formulation. Some extreme cases of ID unsteady shock tube 
and steady nozzle problems are presented in this paper; the performance of these split 
fluxes are compared against the exact solution. The result of accounting for real-gas effect 
is discussed. 

The Euler equations are integrated using the explicit Lax-Wendroff scheme. To ob- 
tain a crisp and monotone shock representation, the TVD scheme based on the above 
split formulas, as described in [11], is employed, along with the super-Bee limiter [17] for 
steepening of the contact discontinuity. 

The computation domain consists of 200 and 100 equally-spaced intervals for the 
shock-tube and nozzle problems, respectively. 

§4.1 Shock-Tube Problem 

The initial conditions are those used in [8]: 

For 0 < x < 5, 

P 4 = 100 atm, 

T 4 = 9000 K, 
u 4 = 0; 

And 5 < x < 10, 

Pi = 1 atm, 

Ti = 300 K, 

tti = 0. 
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The EOS is generated using the widely referred program by Gordon and McBride [18] 
for equilibrium air in the range 250 K < T < 12000 K, 0.1 atm < p < 100 atm, in which 
17 species are included. With over 3600 sets of state points, a least-square fit for pressure 
is obtained with 20 basis functions of (p,e). The resulting standard deviation does not 
exceed 0.2%, the maximum value occurs, as expected, on the boundary. 

The procedure for finding the “exact” solution basically is similar to that of [7] except 
for differences in the details of the numerical steps. We iterate the system via Newton’s 
procedure; the pressure behind the shock P 2 is iterated until the velocity integrated through 
the rarafaction fan and the velocity obtained from the jump relation across the shock are 
equal at the contact discontinuity. The procedure converges in only a few iterations. 

Figures 1-3 show the numerical results of the Roe, Van Leer and Steger- Warming 
splittings for, respectively, p/p 4 ,tt/u 4 ,p/p 4 and e/e 4 , together with exact solutions for 
real and ideal gases. The jump across the contact is rather large, about one order of 
magnitude. This is a difficult case to calculate, as the initial temperatures differ by a factor 
30 and consequently the compositions of the air are completely different. It has been our 
experience that the TVD scheme can handle large differences in pressure very well, but 
not as well if there also are large temperature differences. Nevertheless, our numerical 
results generally show excellent agreement with the exact solution. A crisp shock profile is 
seen, while the effect of dissipation at the contact discontinuity remains as in the case of 
the ideal gas, but accompanied by slight oscillations. Among all three splittings, the Roe 
scheme seems to give the best results, especially near the contact discontinuity. 

The species molar fractions obtained by the Roe splitting are shown in Fig. 4. The 
sharp peak in the NO molar fraction is a numerical result produced by the smearing of 
internal energy at the contact discontinuity. This is because the NO molar fraction is a 
very strong, non-monotone function of temperature (internal energy), and is most stable 
at some intermediate temperature across the “numerical” contact discontinuity. It is seen 
that a large amount of recovery of O 2 and N 2 from O and N occurs across the contact due 
to the large temperature drop, while there is only a minute change across the shock. This 
is also the case for electrons, N + and NO - . 

The results of excluding the real-gas effect are also displayed for air with the same 
initial pressure and temperature as well as at the same final t (=0.0015 sec). This amounts 
to taking p = 0 Ape for the EOS. Significant dissociation and ionization take place at high 
temperatures in the real gas, thereby yielding a much larger specific internal energy and 
lower density than in the ideal gas. The shock and rarefaction wave in the ideal gas 
are travelling at a slower speed; the shock strength and the jump across the contact 
discontinuity also appear weaker. 

Similar numerical results were obtained by Montagne, Yee and Vinokur [19], based 
on the split-flux formulas of [9]. The parameters of the shock-tube problems solved in [19] 
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are different from those of Figures 1-4 and do not go beyond a temperature jump of a 
factor 10 at the contact surface. Disregarding this difference, the results appear to be of 
comparable accuracy. 

§4.2 Steady Nozzle Problem 

Calculations for steady real-gas flows in ID convergent-divergent and divergent nozzles 
are given along with the exact ideal-gas solutions. The area distribution of the nozzle are 
listed below. 

Convergent-divergent nozzle: 

A{x) = 5.5 -4.5 cos(^^Tr), 4 < x < 10; 

x — 4 

A(x) = 1.2 — 0.2 cos( tt), x < 4. 

4 

Divergent nozzle: 

A{x) = 5.5 + 4.5 tanh(0.7x - 3.5), 0 < x < 10. 

Figures 5-7 show the results p/ptooiP/Pooi “/<*«>> and e/«oo for the convergent-divergent 
nozzle. Excellent agreement is achieved by all three schemes with monotone and sharp 
resolution across the shock, although Steger- Warming’s is slightly more dissipative and 
smears over about 2 interior cells. The kink at the throat resulting from the discontinuity 
of area curvature, A"(x), is resolved well. The ideal-gas shock wave in this case is slightly 
further upstream and the jumps in velocity and energy become weaker. As in the shock- 
tube problem, the high inflow temperature produces significant dissociation and ionization. 
Consequently, the internal energy more than doubles that of the ideal gas. The molar 
fraction of N, O, electrons and NO + in Fig. 8 is decreased by the continued drop in 
temperature associated with the acceleration of the gas up until at the shock, but is 
increased abruptly by the shock. The molar fraction of N 2 , O 2 and NO simply reverses 
the trend due to the conservation of mass. 

Calculations are also carried out for flows in the divergent nozzle; the results are given 
in Figs. 9-12. The shock wave is stronger than the above case, but the flow variations are 
generally similar. In this particular set of area geometry and flow parameters, the real-gas 
and ideal-gas shock locations are essentially identical. 


CONCLUDING REMARKS 

Care in deriving flux-vector and flux-difference splittings for a real gas pays off, as 
evidenced by the quality of the numerical results presented above. This is gratifying, 
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because most formulas derived above, viz. the extensions of Van Leer’s differentiable 
flux-vector splitting and Roe’s flux-difference splitting are close to being mathematically 
unique. The Steger- Warming splitting may be derived differently but we saw neither the 
way nor the need to do so. 

From the numerical results for the shock-tube problem it is clear that the constraint 
on the width of a contact discontinuity is even more stringent for a real gas than for 
an ideal gas. The fictitious discrete states in a numerical profile representing a contact 
discontinuity trigger non-trivial chemistry, disturbing the pressure equilibrium that is so 
easily found for the ideal gas. 
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Figure 1. Shock tube problem, Roe flux-difference splitting. 
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Figure 2. Shock tube problem, Steger-Waxming flux-vector splitting. 
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Figure 3. 'Shock tube problem, Van Leer flux- vector splitting. 
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Figure 4. Shock tube problem, molar fractions of species 
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Figure 5. Convergent-divergent nozzle problem, Roe flux-difference splitting. 
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Figure 6. Convergent-divergent nozzle problem, Steger- Warming flux- vector splitting. 
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Figure 7. Convergent-divergent nozzle problem, Van Leer flux-vector splitting. 
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Figure 8. Convergent-divergent nozzle problem, molar fractions of species. Roe splitting. 
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Figure 10. Divergent nozzle problem, Steger- Warming flux-vector splitting. 
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Figure 11. Divergent nozzle problem. Van Leer flux-vector splitting. 
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